#######################################
# Traffic Results - Crossing Borders Project
#######################################

rm(list=ls())
gc()
library(lfe)
library(data.table)
library(haven)
library(ggplot2)
library(viridis)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(extrafont)
library(here)
library(ggthemr)
library(texreg)
ggthemr("light", layout = "plain")

library(cowplot)

source(here("code", "functions.R"))



# Load and prepare data ---------------------------------------------------
traffic <- read_dta(here("data", "traffic", "zst_balanced_yearly_inputR_all.dta"))
setDT(traffic)

traffic$TR_three = NULL
traffic[ treatment == 1 ,TR_three := "0-15 Minutes"]
traffic[ treatment == 0, TR_three := "15-30 Minutes"]


# Main traffic results ----------------------------------------------------
## Event study -------------------------------------------------------------

for (y in 1997:2015) {
  traffic[  ,paste0("border_",y) := 0]
  traffic[ treatment == 1 & year==y ,paste0("border_",y) := 1]
}

traffic[, "trend" := year-1997]
traffic[, "bordertrend" := treatment*trend]

#This includes linear trends to parallize pre-treatment
mod_traffic_es <- felm(trafficperhouroperative ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 + bordertrend  | as.factor(year) + zst | 0 | zst, data = traffic)

#Without linear trends
mod_traffic_es_notrends = felm(trafficperhouroperative ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015  | as.factor(year) + zst | 0 | zst, data = traffic)

event_study_dat <- data.frame(
  y = coef(mod_traffic_es),
  ymin = coef(mod_traffic_es) - qt(.975, df = df.residual(mod_traffic_es))*sqrt(diag(vcov(mod_traffic_es))),
  ymax = coef(mod_traffic_es) + qt(.975, df = df.residual(mod_traffic_es))*sqrt(diag(vcov(mod_traffic_es))),
  Year = as.numeric(gsub("border_", "", names(coefficients(mod_traffic_es))))
)
event_study_dat <- rbind(setDT(event_study_dat), data.table(c(0), c(0), c(0), c(1999)), use.names = F) %>% filter(complete.cases(Year))

plot_traffic_es <- ggplot(data = event_study_dat, aes(x = Year, y = y, ymin = ymin, ymax = ymax, 
                                                      color = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1))),
                                                      shape = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1)))
                                                      )) + 
  geom_vline(xintercept = c(1999, 2004), linetype = 2, color = "black") +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  geom_pointrange() +
  annotate("text", x = 1997.5, y = 55, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 55, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 55, label = "Free Movement", hjust = 0.5) +
  scale_color_manual(values = cbPalette[c(1,2,8)], name = NULL, labels = c("Free Movement", "Transition", "Pre-Reform"), guide = guide_legend(reverse = T)) +
  scale_shape_manual(values = c(16, 17, 18), 
                     name = NULL, 
                     labels = c("Free Movement", "Transition", "Pre-Reform"), 
                     guide = guide_legend(reverse = T)) +
  ylab("Change in traffic per Hour Operative per Counting Station") +
  theme_minimal() + 
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt")) +
  ggtitle("B. Event-Study Estimates") + 
  scale_x_continuous(breaks = seq(1997, 2015, 2))

## DiD ---------------------------------------------------------------------

#Define variables for DiD
traffic[  ,"transBorder15" := 0]
traffic[year>=2000&year<2004&treatment==1  ,"transBorder15" := 1]
traffic[ , "freeBorder15" := 0]
traffic[year>=2004&treatment==1, "freeBorder15" :=1]

mod_traffic_did <- felm(formula = trafficperhouroperative ~ transBorder15 + freeBorder15 + bordertrend | zst + as.factor(year) | 0 | zst, data = traffic)

#Without linear trends
mod_traffic_did_notrends <- felm(formula = trafficperhouroperative ~ transBorder15 + freeBorder15 | zst + as.factor(year) | 0 | zst, data = traffic)

plot_traffic_did <- ggplot(coef.gg.prep(list(mod_traffic_did)), aes(x = modelcoef, y = -1*y, xmin = ylo95, xmax = yhi95, color = variable, shape = variable, group = variable)) +
   geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
   ggstance::geom_linerangeh(position = ggstance::position_dodgev(.25)) +
   ggstance::geom_pointrangeh(aes(xmin = ylo90, xmax = yhi90), size = 1.25, position = ggstance::position_dodgev(.25)) +
   ylab(NULL) +
   scale_y_continuous(breaks = c(-2, -1), labels = c("Transition\nPeriod", "Free\nMovement"), expand = expand_scale(mult = 0, add = 0.5)) +
   scale_color_manual(name = NULL, guide = "none", values = cbPalette[-3]) +
   scale_shape_discrete(name = NULL, guide = "none") +
   xlab("Effect of Border Proximity on Traffic per Operative Hour") +
   ggtitle("C. Difference-in-Differences Estimates") +
   theme_bw() + 
   theme(legend.position = "bottom",
         legend.background = element_rect(color = "black"),
         legend.margin = margin(4,4,4,4,unit = "pt"))


## Time trends -------------------------------------------------------------
plot_traffic_raw <- traffic %>% filter(year >= 1991) %>% 
  ggplot(aes(x = year, y = trafficperhouroperative, color = TR_three, fill = TR_three, linetype = TR_three, shape = TR_three)) + 
  stat_summary(fun.y = mean, geom = "point") + 
  geom_smooth(method = "loess") + 
  geom_vline(xintercept = c(1999, 2004), linetype = "dashed", color = "black") +
  ylab("Traffic per Hour Operative per Counting Station") +
  xlab("Year") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019, 2), expand = expansion(mult = 0, add = 0.5)) +
  scale_color_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_fill_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_linetype_manual(values = c("solid", "dashed", "dotted"), name = NULL, guide = "legend") +
  scale_shape_manual(values = c(16, 17, 18), name = NULL, guide = "legend") +
  scale_y_continuous(breaks = seq(0, 1000, 300)) +
  annotate("text", x = 1997.5, y = 1200, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 1200, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 1200, label = "Free Movement", hjust = 0.5) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"),
        axis.title.y = element_text(vjust = 0)) +
  ggtitle("A. Traffic per Hour Operative in Border Region")



## Output plots
plots <- align_patches(plot_traffic_raw, plot_traffic_did)
top_row <- plot_grid(plots[[1]], plot_traffic_es, align = "h")
pdf(here("figures", "Fig8.pdf"), height = 8.27, width = 12.69)
plot_grid(top_row, plots[[2]], ncol = 1, rel_heights = c(3, 1.5))
dev.off()


# No Nationalstrassen -----------------------------------------------------
traffic_nostr <- read_dta(here("data", "traffic", "zst_balanced_yearly_inputR_noNatStr.dta"))
setDT(traffic_nostr)
traffic_nostr[ treatment == 1 ,TR_three := "0-15 Minutes"]
traffic_nostr[ treatment == 0, TR_three := "15-30 Minutes"]


## Event Study-----------------------------------------------------
for (y in 1997:2015) {
  traffic_nostr[  ,paste0("border_",y) := 0]
  traffic_nostr[ treatment == 1 & year==y ,paste0("border_",y) := 1]
}

traffic_nostr[, "trend" := year-1997]
traffic_nostr[, "bordertrend" := treatment*trend]

mod_traffic_nostr_es = felm(trafficperhouroperative ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 + bordertrend  | as.factor(year) + zst | 0 | zst , data = traffic_nostr)

mod_traffic_nostr_es_notrends = felm(trafficperhouroperative ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015  | as.factor(year) + zst | 0 | zst , data = traffic_nostr)

event_study_dat <- data.frame(
  y = coef(mod_traffic_nostr_es),
  ymin = coef(mod_traffic_nostr_es) - qt(.975, df = df.residual(mod_traffic_nostr_es))*sqrt(diag(vcov(mod_traffic_nostr_es))),
  ymax = coef(mod_traffic_nostr_es) + qt(.975, df = df.residual(mod_traffic_nostr_es))*sqrt(diag(vcov(mod_traffic_nostr_es))),
  Year = as.numeric(gsub("border_", "", names(coefficients(mod_traffic_nostr_es))))
)

event_study_dat <- rbind(setDT(event_study_dat), data.table(c(0), c(0), c(0), c(1999)), use.names = F) %>% filter(complete.cases(Year))

plot_traffic_nostr_es <- ggplot(data = event_study_dat, aes(x = Year, y = y, ymin = ymin, ymax = ymax, 
                                                            color = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1))),
                                                            shape = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1)))
                                                            )) + 
  geom_vline(xintercept = c(1999, 2004), linetype = 2, color = "black") +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  geom_pointrange() +
  annotate("text", x = 1997.5, y = 55, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 55, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 55, label = "Free Movement", hjust = 0.5) +
  scale_color_manual(values = cbPalette[c(1,2,8)], name = NULL, labels = c("Free Movement", "Transition", "Pre-Reform"), guide = guide_legend(reverse = T)) +
  scale_shape_manual(values = c(16, 17, 18), 
                     name = NULL, 
                     labels = c("Free Movement", "Transition", "Pre-Reform"), 
                     guide = guide_legend(reverse = T)) +
  ylab("Change in traffic per Hour Operative per Counting Station") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019,2)) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt")) +
  ggtitle("B. Event-Study Estimates")



## DiD ---------------------------------------------------------------------
traffic_nostr[  ,"transBorder15" := 0]
traffic_nostr[year>=2000&year<2004&treatment==1  ,"transBorder15" := 1]
traffic_nostr[ , "freeBorder15" := 0]
traffic_nostr[year>=2004&treatment==1, "freeBorder15" :=1]

mod_traffic_nostr_did <- felm(formula = trafficperhouroperative ~ transBorder15 + freeBorder15 + bordertrend | zst + as.factor(year) | 0 | zst, data = traffic_nostr)

mod_traffic_nostr_did_notrends <- felm(formula = trafficperhouroperative ~ transBorder15 + freeBorder15 | zst + as.factor(year) | 0 | zst, data = traffic_nostr)

plot_traffic_nostr_did <- ggplot(coef.gg.prep(list(mod_traffic_nostr_did)), aes(x = modelcoef, y = -1*y, xmin = ylo95, xmax = yhi95, color = variable, shape = variable, group = variable)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  ggstance::geom_linerangeh(position = ggstance::position_dodgev(.25)) +
  ggstance::geom_pointrangeh(aes(xmin = ylo90, xmax = yhi90), size = 1.25, position = ggstance::position_dodgev(.25)) +
  ylab(NULL) +
  scale_y_continuous(breaks = c(-2, -1), labels = c("Transition\nPeriod", "Free\nMovement"), expand = expand_scale(mult = 0, add = 0.5)) +
  scale_color_manual(name = NULL, guide = "none", values = cbPalette[-3]) +
  scale_shape_discrete(name = NULL, guide = "none") +
  xlab("Effect of Border Proximity on Traffic per Operative Hour") +
  ggtitle("C. Difference-in-Differences Estimates") +
  theme_bw() + 
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"))

## Time trends-----------------------------------------------------
plot_traffic_nostr_raw <- traffic_nostr %>% filter(year >= 1991) %>% 
  ggplot(aes(x = year, y = trafficperhouroperative, color = TR_three, fill = TR_three, linetype = TR_three, shape = TR_three)) + 
  stat_summary(fun.y = mean, geom = "point") + 
  geom_smooth(method = "loess") + 
  geom_vline(xintercept = c(1999, 2004), linetype = "dashed", color = "black") +
  ylab("Traffic per Hour Operative per Counting Station") +
  xlab("Year") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019, 2), expand = expansion(mult = 0, add = 0.5)) +
  scale_color_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_fill_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_linetype_manual(values = c("solid", "dashed", "dotted"), name = NULL, guide = "legend") +
  scale_shape_manual(values = c(16, 17, 18), name = NULL, guide = "legend") + 
  scale_y_continuous(breaks = seq(0, 1000, 300)) +
  annotate("text", x = 1997.5, y = 1200, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 1200, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 1200, label = "Free Movement", hjust = 0.5) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"),
        axis.title.y = element_text(vjust = 0)) +
  ggtitle("A. Traffic per Hour Operative in Border Region")


plots <- align_patches(plot_traffic_nostr_raw, plot_traffic_nostr_did)
top_row <- plot_grid(plots[[1]], plot_traffic_nostr_es, align = "h")
pdf(here("figures", "Traffic_noNatStr.pdf"), height = 8.27, width = 12.69)
plot_grid(top_row, plots[[2]], ncol = 1, rel_heights = c(3, 1.5))
dev.off()


# Traffic Congestion ------------------------------------------------------
traffic_congestion <- read_dta(here("data", "traffic", "zst_inputR_hourly_conganalysis.dta"))
setDT(traffic_congestion)
traffic_congestion$TR_three = NULL
traffic_congestion[ treatment == 1 ,TR_three := "0-15 Minutes"]
traffic_congestion[ treatment == 0, TR_three := "15-30 Minutes"]

## Event study ------------------------------------------------------
for (y in 1997:2015) {
  traffic_congestion[  ,paste0("border_",y) := 0]
  traffic_congestion[ treatment == 1 & year==y ,paste0("border_",y) := 1]
}

traffic_congestion[, "trend" := year-1997]
traffic_congestion[, "bordertrend" := treatment*trend]

mod_congestion_es <- felm(traffic_cong ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 + bordertrend | as.factor(year) + zst | 0 | zst , data = traffic_congestion)

mod_congestion_es_notrends <- felm(traffic_cong ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 | as.factor(year) + zst | 0 | zst , data = traffic_congestion)

event_study_dat <- data.frame(
  y = coef(mod_congestion_es),
  ymin = coef(mod_congestion_es) - qt(.975, df = df.residual(mod_congestion_es))*sqrt(diag(vcov(mod_congestion_es))),
  ymax = coef(mod_congestion_es) + qt(.975, df = df.residual(mod_congestion_es))*sqrt(diag(vcov(mod_congestion_es))),
  Year = as.numeric(gsub("border_", "", names(coefficients(mod_congestion_es))))
)

event_study_dat <- rbind(setDT(event_study_dat), data.table(c(0), c(0), c(0), c(1999)), use.names = F) %>% filter(complete.cases(Year))

plot_congestion_es <- ggplot(data = event_study_dat, aes(x = Year, y = y, ymin = ymin, ymax = ymax, 
                                                         color = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1))),
                                                         shape = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1)))
                                                         )) + 
  geom_vline(xintercept = c(1999, 2004), linetype = 2, color = "black") +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  geom_pointrange() +
  annotate("text", x = 1997.5, y = 0.05, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 0.05, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 0.05, label = "Free Movement", hjust = 0.5) +
  scale_color_manual(values = cbPalette[c(1,2,8)], name = NULL, labels = c("Free Movement", "Transition", "Pre-Reform"), guide = guide_legend(reverse = T)) +
  scale_shape_manual(values = c(16, 17, 18), 
                     name = NULL, 
                     labels = c("Free Movement", "Transition", "Pre-Reform"), 
                     guide = guide_legend(reverse = T)) + 
  ylab("Change in traffic per Hour Operative per Counting Station") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019,2)) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt")) +
  ggtitle("B. Event-Study Estimates")


## DiD ------------------------------------------------------
traffic_congestion[  ,"transBorder15" := 0]
traffic_congestion[year>=2000&year<2004&treatment==1  ,"transBorder15" := 1]
traffic_congestion[ , "freeBorder15" := 0]
traffic_congestion[year>=2004&treatment==1, "freeBorder15" :=1]

mod_congestion_did <- felm(formula = traffic_cong ~ transBorder15 + freeBorder15 + bordertrend | zst + as.factor(year) | 0 | zst, data = traffic_congestion)
mod_congestion_did_notrends <- felm(formula = traffic_cong ~ transBorder15 + freeBorder15 | zst + as.factor(year) | 0 | zst, data = traffic_congestion)

plot_congestion_did <- ggplot(coef.gg.prep(list(mod_congestion_did)), aes(x = modelcoef, y = -1*y, xmin = ylo95, xmax = yhi95, color = variable, shape = variable, group = variable)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  ggstance::geom_linerangeh(position = ggstance::position_dodgev(.25)) +
  ggstance::geom_pointrangeh(aes(xmin = ylo90, xmax = yhi90), size = 1.25, position = ggstance::position_dodgev(.25)) +
  ylab(NULL) +
  scale_y_continuous(breaks = c(-2, -1), labels = c("Transition\nPeriod", "Free\nMovement"), expand = expand_scale(mult = 0, add = 0.5)) +
  scale_color_manual(name = NULL, guide = "none", values = cbPalette[-3]) +
  scale_shape_discrete(name = NULL, guide = "none") +
  xlab("Effect of Border Proximity on Congestion Index") +
  ggtitle("C. Difference-in-Differences Estimates") +
  theme_bw() + 
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"))

## Time trends ------------------------------------------------------
plot_congestion_raw <- traffic_congestion %>% filter(year >= 1991) %>% 
  ggplot(aes(x = year, y = traffic_cong, color = TR_three, fill = TR_three, linetype = TR_three, shape = TR_three)) + 
  stat_summary(fun.y = mean, geom = "point") + 
  geom_smooth(method = "loess") + 
  geom_vline(xintercept = c(1999, 2004), linetype = "dashed", color = "black") +
  ylab("Congestion Index per Road") +
  xlab("Year") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019, 2), expand = expansion(mult = 0, add = 0.5)) +
  scale_color_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_fill_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_linetype_manual(values = c("solid", "dashed", "dotted"), name = NULL, guide = "legend") +
  scale_shape_manual(values = c(16, 17, 18), name = NULL, guide = "legend") + 
  scale_y_continuous(breaks = seq(0.2, 0.4, 0.05)) +
  annotate("text", x = 1997.5, y = 0.4, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 0.4, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 0.4, label = "Free Movement", hjust = 0.5) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"),
        axis.title.y = element_text(vjust = 0)) +
  ggtitle("A. Congestion Index in Border Region")

## Output results
plots = align_patches(plot_congestion_raw, plot_congestion_did)
top_row = plot_grid(plots[[1]], plot_congestion_es, align = "h")
pdf(here("figures", "TrafficCongestion.pdf"), height = 8.27, width = 12.69)
plot_grid(top_row, plots[[2]], ncol = 1, rel_heights = c(3, 1.5))
dev.off()



# Congestion, No National Strassen --------------------------------------------------------------------
traffic_congestion_nostr <- read_dta(here("data", "traffic", "zst_inputR_hourly_conganalysis_noNatStr.dta"))
setDT(traffic_congestion_nostr)
traffic_congestion_nostr$TR_three = NULL
traffic_congestion_nostr[ treatment == 1 ,TR_three := "0-15 Minutes"]
traffic_congestion_nostr[ treatment == 0, TR_three := "15-30 Minutes"]

## Event study --------------------------------------------------------------------
for (y in 1997:2015) {
  traffic_congestion_nostr[  ,paste0("border_",y) := 0]
  traffic_congestion_nostr[ treatment == 1 & year==y ,paste0("border_",y) := 1]
}

traffic_congestion_nostr[, "trend" := year-1997]
traffic_congestion_nostr[, "bordertrend" := treatment*trend]

mod_congestion_nostr_es <- felm(traffic_cong ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 + bordertrend | as.factor(year) + zst | 0 | zst , data = traffic_congestion_nostr)

mod_congestion_nostr_es_notrends <- felm(traffic_cong ~ border_1997 + 
                   border_1998 + border_2000 + border_2001 + border_2002 + border_2003 + border_2004 + border_2005 + 
                   border_2006 + border_2007 + border_2008 + border_2009 + border_2010 + border_2011 + border_2012 + 
                   border_2013 + border_2014 + border_2015 | as.factor(year) + zst | 0 | zst , data = traffic_congestion_nostr)

event_study_dat <- data.frame(
  y = coef(mod_congestion_nostr_es),
  ymin = coef(mod_congestion_nostr_es) - qt(.975, df = df.residual(mod_congestion_nostr_es))*sqrt(diag(vcov(mod_congestion_nostr_es))),
  ymax = coef(mod_congestion_nostr_es) + qt(.975, df = df.residual(mod_congestion_nostr_es))*sqrt(diag(vcov(mod_congestion_nostr_es))),
  Year = as.numeric(gsub("border_", "", names(coefficients(mod_congestion_nostr_es))))
)

event_study_dat <- rbind(setDT(event_study_dat), data.table(c(0), c(0), c(0), c(1999)), use.names = F) %>% filter(complete.cases(Year))

plot_congestion_nostr_es <- ggplot(data = event_study_dat, aes(x = Year, y = y, ymin = ymin, ymax = ymax, 
                                                               color = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1))),
                                                               shape = factor(ifelse(Year <= 1999, 3, ifelse(Year < 2004, 2, 1)))
                                                               )) + 
  geom_vline(xintercept = c(1999, 2004), linetype = 2, color = "black") +
  geom_hline(yintercept = 0, linetype = 2, color = "black") +
  geom_pointrange() +
  annotate("text", x = 1997.5, y = 0.05, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 0.05, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 0.05, label = "Free Movement", hjust = 0.5) +
  scale_color_manual(values = cbPalette[c(1,2,8)], name = NULL, labels = c("Free Movement", "Transition", "Pre-Reform"), guide = guide_legend(reverse = T)) +
  scale_shape_manual(values = c(16, 17, 18), 
                     name = NULL, 
                     labels = c("Free Movement", "Transition", "Pre-Reform"), 
                     guide = guide_legend(reverse = T)) + 
  ylab("Change in traffic per hour operative per Counting Station") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019,2)) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt")) +
  ggtitle("B. Event-Study Estimates")


## DiD --------------------------------------------------------------------
traffic_congestion_nostr[  ,"transBorder15" := 0]
traffic_congestion_nostr[year>=2000&year<2004&treatment==1  ,"transBorder15" := 1]
traffic_congestion_nostr[ , "freeBorder15" := 0]
traffic_congestion_nostr[year>=2004&treatment==1, "freeBorder15" :=1]

mod_congestion_nostr_did <- felm(formula = traffic_cong ~ transBorder15 + freeBorder15 + bordertrend | zst + as.factor(year) | 0 | zst, data = traffic_congestion_nostr)
mod_congestion_nostr_did_notrends <- felm(formula = traffic_cong ~ transBorder15 + freeBorder15 | zst + as.factor(year) | 0 | zst, data = traffic_congestion_nostr)

plot_congestion_nostr_did <- ggplot(coef.gg.prep(list(mod_congestion_nostr_did)), aes(x = modelcoef, y = -1*y, xmin = ylo95, xmax = yhi95, color = variable, shape = variable, group = variable)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  ggstance::geom_linerangeh(position = ggstance::position_dodgev(.25)) +
  ggstance::geom_pointrangeh(aes(xmin = ylo90, xmax = yhi90), size = 1.25, position = ggstance::position_dodgev(.25)) +
  ylab(NULL) +
  scale_y_continuous(breaks = c(-2, -1), labels = c("Transition\nPeriod", "Free\nMovement"), expand = expand_scale(mult = 0, add = 0.5)) +
  scale_color_manual(name = NULL, guide = "none", values = cbPalette[-3]) +
  scale_shape_discrete(name = NULL, guide = "none") +
  xlab("Effect of Border Proximity on Congestion Index") +
  ggtitle("C. Difference-in-Differences Estimates") +
  theme_bw() + 
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"))

## Time trends --------------------------------------------------------------------
plot_congestion_nostr_raw <- traffic_congestion_nostr %>% filter(year >= 1991) %>% 
  ggplot(aes(x = year, y = traffic_cong, color = TR_three, fill = TR_three, linetype = TR_three, shape = TR_three)) + 
  stat_summary(fun.y = mean, geom = "point") + 
  geom_smooth(method = "loess") + 
  geom_vline(xintercept = c(1999, 2004), linetype = "dashed", color = "black") +
  ylab("Congestion Index per Road") +
  xlab("Year") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(1991, 2019, 2), expand = expansion(mult = 0, add = 0.5)) +
  scale_color_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_fill_manual(values = cbPalette[-3], name = NULL, guide = "legend") +
  scale_linetype_manual(values = c("solid", "dashed", "dotted"), name = NULL, guide = "legend") +
  scale_shape_manual(values = c(16, 17, 18), name = NULL, guide = "legend") +
  scale_y_continuous(breaks = seq(0.2, 0.4, 0.05)) +
  annotate("text", x = 1997.5, y = 0.475, label = "Pre-Reform", hjust = 0.5) +
  annotate("text", x = 2001.5, y = 0.475, label = "Transition", hjust = 0.5) +
  annotate("text", x = 2009.75, y = 0.475, label = "Free Movement", hjust = 0.5) +
  theme(legend.position = "bottom",
        legend.background = element_rect(color = "black"),
        legend.margin = margin(4,4,4,4,unit = "pt"),
        axis.title.y = element_text(vjust = 0)) +
  ggtitle("A. Congestion Index in Border Region")

plots <- align_patches(plot_congestion_nostr_raw, plot_congestion_nostr_did)
top_row <- plot_grid(plots[[1]], plot_congestion_nostr_es, align = "h")
pdf(here("figures", "Congestion_noNatStr.pdf"), height = 8.27, width = 12.69)
plot_grid(top_row, plots[[2]], ncol = 1, rel_heights = c(3, 1.5))
dev.off()


# Regression Tables -------------------------------------------------------

## Traffic results
Texreg(
  list(
    mod_traffic_es,
    mod_traffic_es_notrends,
    mod_traffic_did,
    mod_traffic_did_notrends,
    mod_traffic_nostr_es,
    mod_traffic_nostr_es_notrends,
    mod_traffic_nostr_did,
    mod_traffic_nostr_did_notrends
  ),
  custom.header = list("Traffic per Operative Hour" = 1:8),
  custom.gof.names = c("Num. obs.", "Num.groups: year", "Num.groups: traffic station"),
  custom.gof.rows = list("Linear Trend" = rep(c("Yes", "No"), 4),
                         "Highways Included" = c(rep("Yes", 4), rep("No", 4))),
  custom.coef.map = list(
    "transBorder15" = "0--15 Minutes X Transition",
    "freeBorder15" = "0--15 Minutes X Free",
    "border_1997" = "0--15 Minutes X 1997",
    "border_1998" = "0--15 Minutes X 1998",
    "border_2000" = "0--15 Minutes X 2000",
    "border_2001" = "0--15 Minutes X 2001",
    "border_2002" = "0--15 Minutes X 2002",
    "border_2003" = "0--15 Minutes X 2003",
    "border_2004" = "0--15 Minutes X 2004",
    "border_2005" = "0--15 Minutes X 2005",
    "border_2006" = "0--15 Minutes X 2006",
    "border_2007" = "0--15 Minutes X 2007",
    "border_2008" = "0--15 Minutes X 2008",
    "border_2009" = "0--15 Minutes X 2009",
    "border_2010" = "0--15 Minutes X 2010",
    "border_2011" = "0--15 Minutes X 2011",
    "border_2012" = "0--15 Minutes X 2012",
    "border_2013" = "0--15 Minutes X 2013",
    "border_2014" = "0--15 Minutes X 2014",
    "border_2015" = "0--15 Minutes X 2015",
    "bordertrend" = "0--15 Minutes X (year-1997)"
  ),
  caption = "Effects on Traffic. Outcome: Traffic per Operative Hour. Sample: Border Region, Travel Minutes $<=$ 30, 1997--2015. Baseline Year: 1999. Note: All regressions include year and traffic-station fixed effects.",
  label = "tab:traffic_operative_hour",
  file = here("tables", "traffic_operative_hour.tex"),
  single.row = T, 
  sideways = F,
  fontsize = "scriptsize"
)





## Congestion results
Texreg(
  list(
    mod_congestion_es,
    mod_congestion_es_notrends,
    mod_congestion_did,
    mod_congestion_did_notrends,
    mod_congestion_nostr_es,
    mod_congestion_nostr_es_notrends,
    mod_congestion_nostr_did,
    mod_congestion_nostr_did_notrends
  ),
  custom.header = list("Traffic Congestion" = 1:8),
  custom.gof.names = c("Num. obs.", "Num.groups: year", "Num.groups: traffic station"),
  custom.gof.rows = list("Linear Trend" = rep(c("Yes", "No"), 4),
                         "Highways Included" = c(rep("Yes", 4), rep("No", 4))),
  custom.coef.map = list(
    "transBorder15" = "0--15 Minutes X Transition",
    "freeBorder15" = "0--15 Minutes X Free",
    "border_1997" = "0--15 Minutes X 1997",
    "border_1998" = "0--15 Minutes X 1998",
    "border_2000" = "0--15 Minutes X 2000",
    "border_2001" = "0--15 Minutes X 2001",
    "border_2002" = "0--15 Minutes X 2002",
    "border_2003" = "0--15 Minutes X 2003",
    "border_2004" = "0--15 Minutes X 2004",
    "border_2005" = "0--15 Minutes X 2005",
    "border_2006" = "0--15 Minutes X 2006",
    "border_2007" = "0--15 Minutes X 2007",
    "border_2008" = "0--15 Minutes X 2008",
    "border_2009" = "0--15 Minutes X 2009",
    "border_2010" = "0--15 Minutes X 2010",
    "border_2011" = "0--15 Minutes X 2011",
    "border_2012" = "0--15 Minutes X 2012",
    "border_2013" = "0--15 Minutes X 2013",
    "border_2014" = "0--15 Minutes X 2014",
    "border_2015" = "0--15 Minutes X 2015",
    "bordertrend" = "0--15 Minutes X (year-1997)"
  ),
  caption = "Effects on Traffic Congestion. Outcome: Congestion Index, defined as the absolute difference in traffic between directions of a counting station per hour, divided by the station's mean daily traffic. Sample: Border Region, Travel Minutes $<=$ 30, 1997--2015. Baseline Year: 1999. Note: All regressions include year and traffic-station fixed effects.",
  label = "tab:traffic_congestion",
  file = here("tables", "traffic_congestion.tex"),
  single.row = T, 
  sideways = F, 
  fontsize = "scriptsize"
)

